
rm(list=ls())
library(foreign)
library(survival)

tdata <- data.frame(rx=1)

## First comparison, basic elapsed 

#def.par <- par(no.readonly = TRUE) # save default, for resetting
pdf("test1.pdf",width=9,height=6)
par(mfrow=c(1,2))
#layout(matrix(c(1,1,2,2),2,1, byrow = TRUE))
#par(mar=c(4,5,2,2)+.1, cex.main=1)
sim1 <- read.dta("elapsed12.dta")
fit1 <-coxph(formula = Surv(gapstart,gapstop,status,type = "counting")~ rx + strata(strata), data=sim1, na.action = na.exclude, eps = 0.0001, iter.max = 100, method = "efron")
summary(fit1)
cox.zph(fit1)

surv1 <- survfit(fit1,newdata=tdata)
plot(surv1,fun='cumhaz', lty=1,xlab="Time Since Last Event", ylab="Cumulative Hazard",main="Cumulative Hazard Plot")
betahat <- coef(fit1)
#text(6,1,expression(paste(hat(beta)," = -.867 (.070)")))
#par(mar=c(4,2,2,2)+.1, cex.main=1)
fit1.coxsnell <- abs(sim1$status-fit1$residuals)
surv.fit1.coxsnell <- survfit(Surv(fit1.coxsnell,sim1$status)~1)
plot(surv.fit1.coxsnell$time, -log(surv.fit1.coxsnell$surv),main="Cox-Snell Residual Plot", xlab="Cox-Snell Residuals", ylab="Cumulative Hazard of Cox-Snell Residuals",type="S")
rug(surv.fit1.coxsnell$time)
abline(a=0,b=1)
dev.off()
#par(def.par)
